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LOCAL DISCONTINUOUS GALERKIN METHODS FOR PARTIAL DIFFERENTIAL 
EQUATIONS WITH HIGHER ORDER DERIVATIVES* 

JUE YAN* AND CHI-WANG SHU* 

Abstract. In this paper we review the existing and develop new local discontinuous Galerkin methods 
for solving time dependent partial differential equations with higher order derivatives in one and multiple 
space dimensions. We review local discontinuous Galerkin methods for convection diffusion equations in- 
volving second derivatives and for KdV type equations involving third derivatives. We then develop new 
local discontinuous Galerkin methods for the time dependent bi-harmonic type equations involving fourth 
derivatives, and partial differential equations involving fifth derivatives. For these new methods we present 
correct interface numerical fluxes and prove L 2 stability for general nonlinear problems. Preliminary numer- 
ical examples are shown to illustrate these methods. Finally, we present new results on a post-processing 
technique, originally designed for methods with good negative-order error estimates, on the local discontin- 
uous Galerkin methods applied to equations with higher derivatives. Numerical experiments show that this 
technique works as well for the new higher derivative cases, in effectively doubling the rate of convergence 
with negligible additional computational cost, for linear as well as some nonlinear problems, with a local 
uniform mesh. 

Key words, discontinuous Galerkin method, partial differential equations with higher derivatives, 
stability, error estimate, post-processing 

Subject classification. Applied and Numerical Mathematics 

1. Introduction. In this paper we review existing and develop new local discontinuous Galerkin meth- 
ods for solving time dependent partial differential equations with higher order derivatives in one and multiple 
space dimensions. We consider a sequence of such partial differential equations with increasingly higher order 
derivatives. A hyperbolic conservation law 

d 

Ut + ^2 fi{U) Xi - 0 O- 1 ) 

<= 1 

is a partial differential equation with first derivatives. The convection diffusion equation 

Ut + £ fi(U) Xi - £ = 0 (1.2) 

t= 1 *=1 j = 1 

where ( aij(U )) is a symmetric, semi-positive definite matrix, is a partial differential equation with second 
derivatives. The general KdV type equation 

Ut + ^ fi(U) Xi + ^ | r i(U) ^ 9ij{ r i{U)xi)x j J — 0 (1-3) 

i=l \ 3=1 / x . 
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is a partial differential equation with third derivatives. The time dependent bi-harmonic equation 

d d 

kf “b y! fi{U)xi + )t yr x^x l “0 (1*4) 

2=1 1=1 

is a partial differential equation with fourth derivatives, where the nonlinearity could be more general but 
we just present (1.4) as an example. The following equation 

d d 

L't + 'y " fi(U)xj 4" 9i{UxjX x )xjXjXj = 0 (1.5) 

i~\ i=l 

is a partial differential equation with fifth derivatives, where again the nonlinearity could be more general 
but we just present (1.5) as an example. Similar equations with sixth or higher derivatives could also 
be presented. All these equations, and their time independent steady state counterparts, appear often in 
physical and engineering applications. In this paper we use capital letters U etc. to denote the solutions to 
the PDEs and lower case letters to denote the numerical solutions. 

The type of discontinuous Galerkin methods we will discuss in this paper, using a discontinuous Galerkin 
finite element approximation in the spatial variables coupled with explicit, nonlinearly stable high order 
Runge-Kutta time discretization [19], were first developed for the conservation laws (1.1) containing first 
derivatives by Cockburn et al. in a series of papers [8, 9, 6, 4, 10]. We will briefly review this method in 
section 2.1. For a detailed description of the method as well as its implementation and applications, we 
refer the readers to the lecture notes [3], the survey paper [5], other papers in that Springer volume, and the 
review paper [12]. 

For equations containing higher order spatial derivatives, discontinuous Galerkin methods cannot be 
directly applied. This is because the solution space, which consists of piecewise polynomials discontinuous at 
the element interfaces, is not regular enough to handle higher derivatives. This is a typical “non-conforming” 
case in finite elements. A naive and careless application of the discontinuous Galerkin method directly to the 
heat equation containing second derivatives could yield a method which behaves nicely in the computation 
blit is inconsistent with the original equation and has 0(1) errors to the exact solution [18, 12]. 

The idea of local discontinuous Galerkin methods for time dependent partial differential equations with 
higher derivatives is to rewrite the equation into a first order system, then apply the discontinuous Galerkin 
method on the system. A key ingredient for the success of such methods is the correct design of interface 
numerical fluxes. These fluxes must be designed to guarantee stability and local solvability of all the auxiliary 
variables introduced to approximate the derivatives of the solution. The local solvability of all the auxiliary 
variables is why the method is called a “local” discontinuous Galerkin method in [11]. 

The first local discontinuous Galerkin method was developed by Cockburn and Shu [11], for the con- 
vection diffusion equation (1.2) containing second derivatives. Their work was motivated by the successful 
numerical experiments of Bassi and Rebay [1] for the compressible Navier-Stokes equations. Later, Yan and 
Shu [20] developed a local discontinuous Galerkin method for the general KdV type equation (1.3) containing 
third derivatives. In both [11] and [20], suitable numerical fluxes at element interfaces were given, which led 
to provable nonlinear L 2 stability of the methods as well as error estimates for the linear cases. Numerical 
examples were shown to illustrate the stability and accuracy of these methods. These results will be briefly 
reviewed in sections 2.2 and 2.3, to motivate the new development in later sections. 

In section 3 of this paper we develop new local discontinuous Galerkin methods for the time dependent 
bi-harmonic type equation (1.4) involving fourth derivatives, and in section 4 we do the same thing for the 
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partial differential equation (1.5) involving fifth derivatives. Similar methods can be designed for well posed 
partial differential equations involving even higher derivatives. As before, we give recipes for correct inter- 
element numerical fluxes which lead to provable nonlinear L 2 stability of the methods. The stability result 
is independent of the coefficients in front of the higher order derivatives, hence the methods are especially 
suitable for the so-called “convection dominated” problems, i.e. those with small coefficients to the higher 
derivative terms and hence are dominated by the first derivative convection terms. Also, these methods are 
extremely local and hence efficient for parallel implementations and easy for h-p adaptivity. We will provide 
preliminary numerical examples to verify the stability and accuracy of these methods. 

A post-processing technique was introduced in [7], which is a local convolution if the mesh is locally 
uniform and has been proven to be able to recover the solution to an accuracy of Ax 2k + l rather than the 
usual A:r* +1 when piecewise polynomials of degree k are used, for linear conservation laws (1.1) (theoretically 
proven) and for linear convection diffusion equations (1.2) (numerically observed). The post-processing is 
even useful in enhancing accuracies for some nonlinear cases in numerical experiments. The key ingredient 
of this post-processing technique is a negative-order error estimate, thus if the problem is linear, the solution 
is sufficiently smooth, and the method has a good negative-order error estimate, then the post-processing 
will achieve the order of convergence of the negative-order norm (which is usually bigger than the order of 
convergence of the L 2 -norm of the error in Galerkin methods). In section 5 of this paper, we apply this 
post-processing technique to our new local discontinuous Galerkin method in [20] and in this paper, for the 
linear KdV like equations (1.3), linear time dependent bi-harmonic equations (1.4), and the linear PDE (1.5) 
containing fifth derivatives. We observe the same accuracy enhancement capability for all these cases, where 
accuracy has been increased From Ax k+1 before post-processing to Ax 2fe+1 or higher (in many cases Az 2 * +2 ) 
after post-processing. This strongly suggests that the LDG methods we have designed for these equations 
with higher derivatives have an order of convergence of 2k H- 1 or higher in the negative-order norm, when 
piecewise polynomials of degree k are used. The post-processing technique thus is fully taking advantage of 
this. To demonstrate that this accuracy enhancement can even be observed for some nonlinear problems, 
we also apply it to the local discontinuous Galerkin method in [20] for the nonlinear KdV equations. As the 
post-processing step is applied only at the end of the computation and is applied only locally, the additional 
computational cost is negligible. 

Concluding remarks summarizing results in this paper and indicating future work are included in section 

6 . 

2. Review of the discontinuous Galerkin methods for PDEs with first, second and third 
derivatives. In this section we review the essential points of the discontinuous Galerkin methods for the 
conservation laws (1.1) with first derivatives, the local discontinuous Galerkin methods for the convection 
diffusion equations (1.2) with second derivatives and the KdV type equations (1.3) with third derivatives. 
For simplicity of presentation, we will use one dimensional cases to present the methods. However, we will 
indicate which results are also valid for the general multi-dimensional cases. 

2.1. Conservation Laws. The one dimensional version of the conservation laws (1.1) has the form 

U, + f(U) x = 0 . ( 2 - 1 ) 

First let’s introduce some notations. The computational mesh is given by Ij =[Xj_i ,Xj + i], for j — 
1,..., jV, with the center of the cell denoted by xj = \ (xj_ i + x J+ ^ ) and the size of each cell by A Xj = 
x . . _ x . , We will denote Ax = max, A Xj. If we multiply (2.1) by an arbitrary test function V(x), 

J ~ 2 J 
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( 2 . 2 ) 


integrate over the interval Ij, and integrate by parts, we get 


[ V t Vdx- f f(U)V x dx + f(U{Xj + i,t))V(xj + i) — f(U(xj_i,t))V(Xj_i) = 0. 

Jlj Jlj 


This is the starting point for designing the discontinuous Galerkin method. The semi-discrete version of the 
discontinuous Galerkin method [9, 6, 4, 10] can be described as follows: we replace both the solution U and 
the test function V by piecewise polynomials of degree at most k, and denote them by u and v. That is, 
u,v 6 V&x where 


Vai = {v : v is a polynomial of degree at most k for x € Ij, j = 1 , .... A r } . (2.3) 


With this choice, there is an ambiguity in (2.2) in the last two terms involving the boundary values at 
x j± i , as both the solution u and the test function v are discontinuous exactly at these boundary points. A 
crucial ingredient for the success of discontinuous Galerkin method for the conservation laws is the correct 
choice of the numerical fluxes to overcome (or we could say it in a positive way, to utilize) this ambiguity 
of discontinuity at the element interfaces. The idea is to treat these terms by an upwinding mechanism 
(information from characteristics), borrowed from the successful high resolution finite volume schemes. Thus 
f(u) at the interface Xj+i for each j is given by a single valued monotone numerical flux 




i+r 


(2.4) 


which depends both on the left limit u j and on the right limit ut, j of the discontinuous numerical solution 

J-r ? 

at the interface x j+ 1 . Here monotone flux means that the function / is a non-decreasing function of its first 
argument and a non-increasing function of its second argument. It is also assumed to be at least Lipschitz 
continuous with respect to each argument and to be consistent with the physical flux f(u) in the sense that 
f(u,u) = f(u). If (2.1) is a system rather than a scalar equation, then the monotone flux is replaced by an 
exact or approximate Riemann solver, see, e.g. [16] for details. On the other hand, the test function v at the 
interfaces Xj± 1 is taken from inside the cell Ij, namely v~ , and v 1 , respectively. The final semi-discrete 
scheme 


/ utvdx — / 

Ji 3 Ji 3 


f{u)v x dx + 4+ 


h v j+$ 


~ U- 


i v : 


= 0 


(2.5) 


is then discretized by a nonlinearly stable high order Runge-Kutta time discretizations [19]. Nonlinear TVB 
limiters [17] may be used if the solution contains strong discontinuities. The schemes thus obtained have the 
following attractive properties, for the general multi-dimensional case (1.1) with arbitrary triangulations: 

1. It can be easily designed for any order of accuracy. In fact, the order of accuracy can be locally 
determined in each cell, thus allowing for efficient p adaptivity. 

2. It can be used on arbitrary triangulations, even those with hanging nodes, thus allowing for efficient 
h adaptivity. 

3. It is extremely local in data communications. The evolution of the solution in each cell needs to 
communicate only with the immediate neighbors, regardless of the order of accuracy, thus allowing 
for efficient parallel implementations. See, e.g. [2]. 

These schemes also have the following provable theoretical properties, all of these are valid also for the 
multi-dimensional case (1.1) with arbitrary triangulations: 

1. The semi-discrete scheme (2.5), and certain time discretization of it, such as implicit backward 
Euler and Crank-Nicholson, have excellent nonlinear stability properties. One can prove a strong L' 2 
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stability and a cell entropy inequality for the square entropy, for the general nonlinear scalar case 
(1.1), for any orders of accuracy on arbitrary triangulations in any space dimension, without the 
need for nonlinear limiters [14]. Notice that these stability results are valid even when the solution 
contains discontinuities such as shocks. 

2. For linear problems with smooth solutions, these methods using piecewise polynomials of degree k 
have a provable error estimate of order Ax k+ * in L 2 [15]. In effect, for most triangulations one 
could observe (and prove in many cases) convergence of the order Ax k+1 in both L 2 and L°° norms, 
for both linear and nonlinear problems. 

3. When nonlinear TVB limiters [17, 9, 4] are used, the methods can be proven stable in the total 
variation norm for scalar one dimensional nonlinear problems (2.1), and stable in the L norm for 
scalar multi-dimensional nonlinear problems (1.1). 

2.2. Convection diffusion equations. The one dimensional version of the convection diffusion equa- 
tion (1.2) has the form 

U t + f(U) x - (a(U)U x ) x = 0, (2 6) 


where a(U ) > 0. 

The semi-discrete version of the local discontinuous Galerkin method for solving (2.6) [11] approximates 
the following lower order system 

Ut + f(U) x - (b(U)Q) x =0, Q - B{U) X = o, (2.7) 


where b(U) = y/rijJ), and B(U) = f U b(U)dU. We can then formally use the same discontinuous Galerkin 
method for the convection equation to solve (2.7), resulting in the following scheme: find u.q 6 Vax suc h 
that, for all test functions v,w € Vax> 


f u t vdx - [ ( f{u ) - b(u)q)v x dx + (/,•+> - b j+ 1 q j+ i )m + i - (/;-* ~ - 0 

Ji, Jlj 

/ qwdx + / B(u)w x dx — B 

Jlj 


+ d i-l w i-k =0 - 


( 2 . 8 ) 


Again, a crucial ingredient for the method to be stable is the correct choice of the numerical fluxes (the 
“hats”). However, there is no longer a upwinding mechanism or characteristics to guide the design of these 
fluxes. In [11], criteria are given for these fluxes to guarantee stability and convergence. The best choice is 
to use the “alternating principle” in designing the fluxes: 


b = 


B(u+) - B{u 


!/+ _ ' 


B = B{u + ) 


(2.9) 


and / is chosen as before in (2.4). Notice that we did not write the subscript j + \ for the fluxes as they are 
all evaluated at this interface point. The “alternating principle” refers to the alternating choices of q and B: 
if the left value is chosen for the former then the right value is chosen for the latter, as in (2.9). One could 
also choose 


■ + ; B = B(u~) 


with all the other fluxes unchanged. These choices of fluxes guarantee stability and optimal convergence. 

We remark that the appearance of the auxiliary variable q is superficial: when a local basis is chosen 
in cell Ij then q is eliminated by using the second equation in (2.8) and solving a small linear system if the 
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local basis is not orthogonal. The actual scheme for u takes a form similar to that for convection alone. This 
is a big advantage of the scheme over the traditional “mixed methods”, whose auxiliary variable is usually 
genuinely global. This is the reason that the scheme is termed local discontinuous Galerkin method in [11]. 
We also remark that the choice of fluxes in (2.9) by the alternating principle yields the most compact stencil 
for the scheme for u after the auxiliary variable q is locally eliminated. 

The schemes thus designed for the one dimensional equation (2.6), or in fact for the most general multi 
dimensional nonlinear convection diffusion equations (1.2), which is nonlinear both in the first derivative 
convection part and in the second derivation diffusion part, retain all of the three attractive properties 
listed above for the method used on convection equations. They also have the following provable theoretical 
properties, all of these are valid for the multi-dimensional case (1.2) with arbitrary triangulations [11]: 

1. The semi-discrete scheme (2.8), and certain time discretization of it, such as implicit backward 
Euler and Crank-Nicholson, have excellent nonlinear stability properties. One can prove a strong 
L 2 stability for the general nonlinear scalar case (1.2), for any orders of accuracy on arbitrary 
triangulations in any space dimension. This stability is independent of the size of the diffusion 
terms and hence is also valid in the limit when the diffusion coefficient goes to zero. 

2. For linear problems with smooth solutions, these methods using piecewise polynomials of degree k 
have a provable error estimate of order Ax k in L 2 . In effect, for most triangulations one could observe 
(and prove in many cases) convergence of the order Ax k + X in both L 2 and L°° norms, for both linear 
and nonlinear problems. The same error estimate can be obtained also for the auxiliary variable q, 
which approximates the derivative U x , even if q can be locally eliminated in actual calculation. 

2.3. KdV like equations. The one dimensional version of the KdV like equation (1.3) has the form 

U t + f(U) x + {r'(U)g{r(U) x ) x ) x = 0, (2.10) 

where f(U ), r(U) and g(U) are arbitrary functions. 

The semi-discrete version of the local discontinuous Galerkin method for solving (2.10) [20] approximates 
the following lower order system 


U t + (f(U) 4- r\U)P) x — 0, P - g(Q) x = 0, Q - r(U) x = 0. (2.11) 


We can then formally use the same discontinuous Galerkin method for the convection equation to solve 
(2.10), resulting in the following scheme: find u,p,q 6 Va* such that, for all test functions v,w,z £ Va*, 


/ u t vdx ~ / {f(u) + r f (u)p)v x dx+ (f + r , p\ v.. j - (f + r'fi) x = 0, 

Jlj Jlj V ' j+5 V /j-L J~ 2 

pwdx + h g(q)w x dx - g j+ iw~ + , + gj_ i , = 0, 

I qzdx + / r(u) 

Jlj Jlj 


( 2 . 12 ) 


z x dx — 7\ , i 


4 _ 


'i+iVi +r J-h z I-h 


Again, a crucial ingredient for the method to be stable is the correct choice of the numerical fluxes (the 
“hats”). It is found out in [20] that one can take the following simple choice of fluxes to guarantee stability 
and convergence: 


f = f{u ,u + ), 


r(u + ) — r(u ) 


p = p + , 


9 = 9(Q ,9 + ), r = r{u ) 


(2.13) 
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where f(u~,u+) is a monotone flux for /(«), and -g(q~,q + ) is a monotone flux for -g(q). In fact, the 
crucial part is still the “alternating principle” to take p and f from opposite sides. Thus 

f = f(u~,u + ), ? = — l — , p = p ~ , g = g(q~,q + ), f = r(u + ) 

would also work. Again, the appearance of the auxiliary variables p and q is superficial: when a local 
basis is chosen in cell Ij then both of them can be eliminated by using the second and third equations in 
(2.12) and solving two small linear systems if the local basis is not. orthogonal. The actual scheme for u 
takes a form similar to that for convection alone. We also remark that the choice of fluxes in (2.13) by the 
alternating principle yields a compact stencil for the scheme for u after the auxiliary variables p and q are 
locally eliminated. 

The schemes thus designed for the KdV like equation (2.10), or in fact for the most general multi 
dimensional nonlinear KdV like equations (1.3), which is nonlinear in all the derivatives, retain all of the 
three attractive properties listed above for the method used on convection equations. They also have the 
following provable theoretical properties, all of these are valid for the multi-dimensional case (1.3) with 
arbitrary triangulations [20]: 

1. The semi-discrete scheme (2.12), and certain time discretization of it, such as implicit backward 
Euler and Crank-Nicholson, have excellent nonlinear stability properties. One can prove a strong L 2 
stability and a cell entropy inequality for the square entropy, for the general nonlinear scalar case 
(1.3), for any orders of accuracy on arbitrary triangulations in any space dimension, without the 
need for nonlinear limiters [20]. Notice that these stability results are valid even in the limit when 
the coefficients of the dispersive third derivative terms tend to zero. 

2. For one dimensional linear problems with smooth solutions, these methods using piecewise polyno- 
mials of degree k have a provable error estimate of order Ax fc+5 in L 2 [20]. In numerical experiments 
one observes convergence of the order Ax* +1 in both L 2 and T°° norms for both one and multiple 
dimensional linear and nonlinear cases. 

3. A local discontinuous Galerkin method for the bi-harmonic type equations. In this section, 
we present and analyze a LDG method for the bi-harmonic type equation (1.4). We will concentrate on the 
one dimensional case 

Ut + /(to* + (a(U x )U xx ) xx =0, 0 < x < 1 (3.1) 

with an initial condition 

U{x,0) = U°(x), 0 < x < 1 (3-2) 

and periodic boundary conditions. Here a(U x ) > 0. The assumption of periodic boundary conditions is for 
simplicity only and is not essential: the method can be easily designed for non-periodic boundary conditions. 
The generalization to the multiple dimensional case (1.4) is straightforward, following the lines in [11] and 
[20]. 

To define the LDG method, we first introduce the new variables 

R = U X , Q = B(R) X , P = (b(R)Q) x , (3-3) 

where b{R) = y/a{R) and B(R) = J R b{R)dR , and rewrite the equation (3.1) as a first order system: 

U t + if{U) + P) X = o, P - (b(R)Q) x =0, Q - B(R) X =0, R - U x = 0. (3.4) 
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The LDG method is obtained by discretizing the above system with the discontinuous Galerkin method. This 
is achieved by multiplying the four equations in (3.4) by four test functions v,w,z,s respectively, integrate 
over the interval and integrate by parts. We need to pay special attention to the boundary terms resulting 
from the procedure of integration by parts. Thus we seek piecewise polynomial solutions u,p, q, r £ V^x, 
where V^ x is defined in (2.3), such that for all test functions v,w,z,s € V& x we have, for 1 < j < A r , 


Ji utvdx - h ( f{u ) + p)v x dx + (/ + Pj. + , v. + i - (/ + p)^ , i = 
/ pwdx 4- / 

Jij Jij 


pwdx + / b(r)qw x dx — 


/ qzdx -f / B(r)z 


w . i + [b<i ) ! = 0, 


- B 7 , 1 2 ! -fB, iz + ! = 0, 

J+5 3 2 J-k 1 


2 J- 5 


(3.5) 


/ rsdx + / us x dx - u , , i s . , , + fL_ i s + j = 0. 

Ji } ? j -2 


The only ambiguity in the algorithm (3.5) now is the definition of the numerical fluxes (the “hats”), which 
should be designed to ensure stability. It turns out that we can take the simple choices (we omit the subscripts 
j ± \ in the definition of the fluxes as all quantities are evaluated at the interfaces Xj±i) 

/=/(u~iU + ), p = p~, b = ^ i q = <7 + , B = B(r~), u = u+ (3.6) 


where /(u~,u+) is a monotone flux for f(u) in (2.4). Here again, the “alternating principle” in designing 
the fluxes related to the fourth derivatives is at play: the fluxes for p, q , B(r) and u are alternatively taken 
from left and right. We could thus also take 

f = f(u~,u + ), P = p + , b = g(r+) - g(r } , q = q~, B = B(r + ), u = u~ (3.7) 

7 — r 

and stability is still valid. 

We remark again that the appearance of the auxiliary variables p, q and r is superficial: when a local 
basis is chosen in cell Ij then all of them can be eliminated by using the second, third, and fourth equations 
in (3.5) and solving three small linear systems if the local basis is not orthogonal. The actual scheme for u 
takes a form similar to that for convection alone. 

We have the following L 2 stability result for the scheme (3.5)-(3.6). This is similar to the stability result 
in Proposition 2.1 in [11] for the second order convection diffusion equations. 


Proposition 3.1. ( L 2 stability) The solution of the scheme (3.5)-(3.6) satisfies 

li (^-f -) dx+ L ,2<ii<)<<xso 

Proof: We sum up the four equalities in (3.5) and introduce the notation 

Bj(u,p,q,r-,v,w,z,s) = / u t vdx - / (f(u)+p)v x dx+(f + p) v7, , 
Jij Ji , v h+\ 3+ * 

-\f + P) , ut i + / pwdx + / b(r)qw x dx — 

v h-h 1 2 Jij Ji, 

+ j + l, qzdx 4- J B(r)z x dx — B 


tu., , 

i+h J+5 


+Bj_\zt j + / rsdx + / 

3 5 •'~ 5 Ji j Ji, 




rsdx 4- / us x dx — 


(3.8) 


(3.9) 


8 



(3.10) 


Clearly, the solutions u, p, q, r of the scheme (3.5)-(3.6) satisfy 

Bj(u,p,q,r;v,w, z,s) = 0 

for all v,w,z,s £ V'aj- ■ We then take 

v = u, w =r, z = q, s = -p 


to obtain, after some algebraic manipulations, 


0 = Bj(u,p, q, r; u,r,q, —p) 

' u 2 {x,t) 


= Jtf, ( “ ( 2 t ') dX + J, q2i < x ^ dx+ (^i+i +0 i-5 


with the numerical entropy flux H defined by 

H = -F(u~) - p~u~ + B{r~)q~ + (/ + pj u~ - ( bqj r~ - Bq~ + 


up 


(3.11) 


and the extra term 0 given by 

0 = [F( u) +pu- B(r)q] - (/ + p) [u] + bq[r] + B[q] - u[p], 

Here F{u) = /“ f(u)du, and [a] = v + -v~ denotes the jump of a. Notice that we have dropped the subscripts 
about the location j — ± or j + \ as all these quantities are defined at a single interface and depend only on 
the left and right values at that interface. Now all we need to do is to verify 0 > 0. To this end, we notice 
that, with the definition (3.6) of the numerical fluxes and with simple algebraic manipulations, we easily 

obtain 

[pit - B(r)q] - p[u] + bq[r] + B[q ] - u\p] = 0 


and hence 

0 = [F(u)] - f[u] = f (f(s)-f(u~,u + ))ds> 0, 

J U~ 

because of the monotonicity of the flux /. Summing up (3.11) over j would now give the desired L 2 stability 

(3.S). ■ 

For a preliminary numerical example of the methods developed in this section, see Example 5.3 in section 
5. 

4. A local discontinuous Galerkin method for equations with fifth derivatives. In this section, 
we present and analyze a LDG method for the equation (1.5) with fifth derivatives. We will concentrate on 
the one dimensional case 

Ut + f(U) x + g(U X x)xxx — 0 , A x 5 ; 1 ( 4 * 1 ) 


with an initial condition 

U(x,0) = U°(x), 0 < a; < 1 (4-2) 

and periodic boundary conditions. Here g(U X x) is an arbitrary function of U xx . The form of the nonlinearity 
is not general and is chosen simply as an example. The assumption of periodic boundary conditions is also for 
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simplicity only and is not essential: the method can be easily designed for non-periodic boundary conditions. 
The generalization to the multiple dimensional case (1.5) is straightforward, following the lines in [11] and 
[ 20 ]. 

To define the LDG method, we first introduce the new variables 


S = U X , R = S X , Q = g(R) x , P = Q X , 
and rewrite the equation (4.1) as a first order system: 


(4.3) 


Ut + (f(U) + P) x = 0, P-Q x = 0, Q - g(R) x = 0, R-S x = 0, S-U x = 0. (4.4) 


The LDG method is obtained by discretizing the above system with the discontinuous Galerkin method. 
This is achieved by multiplying the five equations in (4.4) by five test functions v,w,y,z,d respectively, 
integrate over the interval I v and integrate by parts. We need to pay special attention to the boundary 
terms resulting from the procedure of integration by parts. Thus we seek piecewise polynomial solutions 
G V Ax , where is defined in (2.3), such that for all test functions v,w,y,z,d € Va? we have, 
for 1 < j < N, 


/ u t vdx- / {f(u)+p)v x dx+ (f + p) u,,i-(7" + p) f + _i=0, 

Jlj Jlj ' ' j+ 5 ■' + 2 V ' j~k 3 2 

J Jm)dx + L qw x dx - q j+ iw j+ ± + ^ = 0, 

<]ydx + g(r)y x dx - g j+ iyj + , +g j _^y+_ l _ =0, 

/ rzdx + / sz x dx - s i+ iz~ , +s. i z + , = 0. 
Jij Ji j 1+2 J+ s J 2 

/ sddx+ / ud x dx - u i+1 d~ t + u,_ i d + ! =0. 

Jlj Jlj 3+2 3 2 J~k 


(4.5) 


The only ambiguity in the algorithm (4.5) now is the definition of the numerical fluxes (the “hats”), which 
should be designed to ensure stability. It turns out that we can again take the simple choices 


/ = /(** > u+ )> p-p ) g = g + , 9 = g(r ,r + ), s = s~, u = u + (4.6) 

where f(u~,u + ) is a monotone flux for f(u) in (2.4), and - 3 (r~,r + ) is a monotone flux for -g{r). Here 
again, the “alternating principle” in designing the fluxes related to the fifth derivatives is at play: the fluxes 
for p, q, s and u are alternatively taken from left and right. We could thus also take 


/ = /(“ » M+ ). p = p + , g = q , g = g(r ,r + ), 


s = s + . 


u = u 


(4.7) 


and stability is still valid. 

As before, the appearance of the auxiliary variables p, q, r and s is superficial: when a local basis is 
chosen in cell Ij then all of them can be eliminated by using the second, third, fourth and fifth equations 
in (4.5) and solving four small linear systems if the local basis is not orthogonal. The actual scheme for u 
takes a form similar to that for convection alone. 

We have the following cell entropy inequality for the square entropy, which implies L 2 stability, when 
summed up over j, for the scheme (4.5)-(4.6). This is similar to the cell entropy inequality in [14] for 
conservation laws and the one in [20] for KdV type equations involving third derivatives. 
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Proposition 4.1. (cell entropy inequality) The solution of the scheme (4.5)-(4.6) satisfies 


-! 

dt J, j 


u 2 (x, t) 


dx 




3 + 5 


(4.8) 


for some numerical entropy flux H J+ i- 

Proof: We sum up the five equalities in (4.5) and introduce the notation 

Bj{u,p,q,r,s;v,w,y,z,d) = J u t vdx - (f(u) +p)v x dx + (/ + p)^ + ^ v jj 

- (/ + p) , vt_ i + Ji pwdx + Ji qw x dx - q j+ iwJ + ± 

» + / 9 V dx + / 9{r)yxdx-g j+ LyJ + i 

2 3 * Jlj Jl 3 

+9j-i»]L i + f' rzdx + l sz * dx ~ S i+h z j+i + *i~i z l-h 

+ / sdtix + / ud x dx - u j+ idj + i 
Jlj Jlj 

Clearly, the solutions u, p, g, r and s of the scheme (4. 5)- (4. 6) satisfy 

Bj (u,p,q,r, S',v,w,y, z,d) = 0 

for all v,w,y,z,d £ Vax- We then take 

v = u, w = s, y = -r, 2 = 9 , d = —p 


(4.9) 


(4.10) 


+ 0 


to obtain, after some algebraic manipulations, 

0 = Bj(u,p, q,r,s:u,s, -r.q, -p) = — ^ ' ) dx + (^+5 _ ^J-i) 

with the numerical entropy flux H defined by 

H = —F(u~) -p~u~ + q~s~ - G(r ~ ) + + p) u~ - 9s - + gr~ - sq~ + up 


J-5 


(4.11) 


and the extra term 0 given by 

0 = [ F(u ) + pu — 9 s + G(r)] — {j + p^j [u] + 9 W — 9[ r ] + — d\p]. 

Here F(u) = f f(u)du and G(r) = f g{r)dr, and [v] = v+ - v~ denotes the jump of v. Now all we need 
to do is to verify 0 > 0. To this end, we notice that, with the definition (4.6) of the numerical fluxes and 
with simple algebraic manipulations, we easily obtain 

\pu - qs ] - p[u\ + 9 [s] + s[q] - u\p\ = 0 


and hence 

0 = [F(u)} - f[u\ - [G(r)] + s[r] 

= f V (/( 0 -/(«"> u+ ))^-^_ {9(0-9(r-,r+))d(, 

> 0 " 

where the last inequality follows from the monotonicity of the fluxes / and —g. This finishes the proof. ■ 
For a preliminary numerical example of the methods developed in this section, see Example 5.4 in section 

5. 
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5. Post-processing and numerical results. In [7], an efficient local post-processing technique has 
been developed, which is applied only at the end of the calculation to the numerical solution at t = T, and 
involves only a local linear operation (convolution) using the information from a fixed (independent of the 
mesh size Ax) number of neighboring cells. The description of this post-processing procedure is omitted 
here because of space limitation. Please see [7] for details. 

It is proven in [7] that this post processing will improve the order of accuracy from Ax k + X to Ax 2 * +1 
when piecewise polynomials of degree k are used, for linear conservation laws (1.1) with a uniform (or locally 
uniform) mesh with smooth solutions. Although not proven theoretically, it is also verified numerically in 
[7] that the same enhancement of accuracy can be observed for linear convection diffusion equations (1.2). 
In fact, the crucial ingredient this post-processing technique depends on is a negative-order error estimate, 
which is usually higher than the order of convergence of the L 2 -norm of the error, for Galerkin type methods. 
If the problem is linear, the solution is sufficiently smooth, then the post-processing will achieve the order 
of convergence of the negative-order norm. 

In this section, we perform numerical experiments of this post-processing techniques when the local 
discontinuous Galerkin method is applied to the linear and nonlinear KdV equations (1.3), and linear versions 
of (1.4) and (1.5) containing fourth and fifth derivatives. We observe the same accuracy enhancement. This 
strongly suggests that the LDG methods we have designed for these equations with higher derivatives have 
an order of convergence of 2/c+ 1 or higher in the negative-order norm, when piecewise polynomials of degree 
k are used. The post-processing technique thus is fully taking advantage of this. Time discretization is by 
the third and fourth order implicit Runge-Kutta methods in [13]. An efficient time marching of these local 
discontinuous Galerkin methods without sacrificing its local property and parallel flexibility is left for future 
work. We have chosen At suitably small so that spatial errors dominate in the numerical results. 

Example 5.1. We solve the linear KdV equation (1.3) in ID: 

Ut + Uxxx = 0 (5.1) 

with initial condition U(: r,0) = sin(ar) and periodic boundary conditions. The exact solution is given by 
U(x,t) = sin(z + t). The numerical errors and order of accuracy before and after post-processing can be 
found in Table 5.1. We clearly observe an accuracy of Ax k + X before the post-processing and Ax 2 * +1 or higher 
(closer to Ax 2 * +2 ) after it, w r hen piecewise P k elements are used. W 7 e do not present the result with k = 0 
because for this case the results before the post-processing is already showm in [20] and the post-processing 
does not increase the order of accuracy (2 A: -f 1 = k + 1 for k = 0). 

Example 5.2. We compute the classical soliton solution of the nonlinear KdV equation (1.3) in ID: 

U t -S(U 2 ) x + U xxx = 0 (5.2) 

with an initial condition u(x,0) = -2sech 2 (x) for -14 < x < 16, see [20] for details. The exact solution 
is given by U{x,t) = — 2sech 2 (x — 4 1). The numerical errors and order of accuracy before and after post- 
processing can be found in Table 5.2. We clearly observe an accuracy of Ax k + l before the post-processing 
and Ax 2k + l or higher after it, wffien piecewise P k elements are used. This indicates that the post-processing 
technique is effective in enhancing orders of accuracy even for this nonlinear PDE. 

Example 5.3. We solve the linear bi-harmonic equation (1.4) in ID: 

Ut T U XX xx ~ 0 (5.3) 
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Table 5.1 

Ut + U x xx = o. U{x, 0) = sin(i). Periodic boundary conditions. L°° errors. Uniform meshes with N cells. LDG methods 
with k = 1,2,3. t - 1. Before (B) and after (A) post-processing. 


k 


N=10 

Z 

II 

to 

o 

N=40 

N=8C 




error 

error 

order 

error 

order 

error 

order 

1 

B 

5.0265E-02 

1.3623E-02 

1.88 

1.0619E-03 

1.96 

8.8570E-04 

1.98 


A 

5.7711E-03 

5.8787E-04 

3.29 

6.3866E-05 

3.20 

7.3564E-06 

3.11 

2 

B 

2.9084E-03 

3.6532E-04 

2.99 

4.6186E-05 

2.98 

5.7816E-06 

2.99 


A 

2.0588E-04 

3.6032E-06 

5.83 

6.3684E-08 

5.82 

1.2763E-09 

5.64 



N=10 

N=20 

N=40 

N=5( 

) 

3 

B 

9.2247E-05 

6.0315E-06 

3.93 

3.8021E-07 

3.98 

1.5583E-07 

3.99 


.4 

2.2816E-05 

9.7519E-08 

7.87 

3.9135E-10 

7.96 

7. 0005 E- 11 

7.71 


Table 5.2 

U t _3 (f/ 2 ) + [/ xxx = 0. I/(x,0) = -2 sech 2 (x). L°° errors. Uniform meshes with N cells. LDG methods with k = 1,2,3. 
t = 0.5. Before (B) and after (A) post-processing. 


k 


o 

ao 

II 

z 

N=160 

N=320 

N=64 

0 



error 

error 

order 

error 

order 

error 

order 

1 

B 

7.1035E-02 

1.4334E-02 

2.30 

4.3454E-03 

1.72 

1.1852E-03 

1.87 


A 

4.8823E-02 

6.0560E-03 

3.01 

7.3566E-04 

3.04 

8.9456E-05 

3.03 

2 

B 

4.8540E-03 

6.2916E-04 

2.94 

8.0172E-05 

2.97 

1.0047E-05 

2.99 


A 

3.0621E-03 

8.0491E-05 

5.24 

1.6184E-06 

5.63 

3.3110E-08^ 

5.61 



o 

00 

II 

z 

N— 160 

N=320 

N=50 

0 

3 

B 

3.5385E-04 

2.4394E-05 

3.85 

1.5343E-06 

3.99 

2.5907E-07 

3.98 


A 

1.6221E-03 

1.6713E-05 

6.60 

8.9620E-08 

7.54 

3.3679E-09 

7.35 


with initial condition U(x, 0) = sin(x) and periodic boundary conditions. The exact solution is given by 
U(x,t) = e _t sin(x). The numerical errors and order of accuracy before and after post-processing can be 
found in Table 5.3. We clearly observe an accuracy of Ax k+i before the post-processing and Ax 2k ' or 
higher (closer to Ax 2k+2 ) after it, when piecewise P k elements are used. We have used multiple precision 
FORTRAN to compute part of this problem, with an effective precision higher than the standard double 
precision, to overcome some of the round-off problems for this problem. If such multiple precision procedure 
is not used, then the post-processed error can only go down to around 10“ 10 level and cannot decrease 
further. The result in Table 5.3 verifies both the stability and accuracy of the local discontinuous Galerkin 
method developed in section 3 and also the effectiveness of the post- processing technique for such methods 
applied to this type of equations. It also strongly suggests that there is a better negative norm error estimate 
on the order of Ax 2fc+2 for piecewise P k elements to this problem. 

Example 5.4. We solve the linear equation (1.5) in ID: 

U t + U xxxxx = 0 (5-4) 

with initial condition C/(x, 0) = sin(x) and periodic boundary conditions. The exact solution is given by 
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Table 5.3 

Ut + Uxxxx — 0. t/(x, 0) = sin(x). Periodic boundary conditions. L°° errors. Uniform meshes with N cells. LOG methods 
with k = 0, 1,2,3. t = 1. Before (B) and after (A) post-processing. 


k 


N=10 

N=20 

o 

■'T 

II 

z 

N=80 



error 

error 

order 

error 

order 

error 

order 

0 

B 

1.1 125E-01 

5.4352E-02 

1.03 

2.7001E-02 

1.00 

1.3478E-02 

1.00 

1 

B 

2.2038E-02 

5.2262E-03 

2.07 

1.3119E-03 

1.99 

3.2831E-04 

1.99 


A 

7.5566E-04 

4.5478E-05 

4.05 

2.8624E-06 

3.98 

1.7930E-07 

3.99 

2 

B 

1.1183E-03 

1.3512E-04 

3.04 

1.6988E-05 

2.99 

2.1265E-06 

2.99 


A 

9.4333E-05 

1.2978E-06 

6.18 

1.9072E-08 

6.08 

2.9261E-10 

6.02 

3 

B 

6.1004E-05 

2.3484E-06 

4.69 

1.4022E-07 

4.06 

8.7476E-09 

4.00 


A 

3.2399E-05 

1.6632E-07 

7.60 

6.9985E-10 

7.89 

2.7864E-12 

7.97 


Table 5.4 

Ut T Uxxxxx = 0. f/(x,0) = sin(x). Periodic boundary conditions . L°° errors. Uniform meshes with N cells. LOG 
methods with k = 0, 1,2,3. t = 1. Before (B) and after (A) post-processing. 


k 


N=10 

N=20 

o 

II 

Z 

z 

II 

Oc 

o 

error 

error 

order 

error 

order 

error 

order 

0 

B 

5.6284E-01 

2.5073E-01 

1.16 

1.1647E-01 

1.10 

5.5926E-02 

1.05 

1 

B 

4.8495E-02 

1.3661E-02 

1.82 

3.5772E-03 

1.93 

8.9887E-04 

1.99 

A 

6.6528E-03 

6.3604E-04 

3.38 

6.2508E-05 

3.34 

6.7327E-06 

3.21 

2 

B 

2.9071E-03 

3.6668E-04 

2.98 

4.6236E-05 

2.98 

5.7831E-06 

2.99 

A 

2.3284E-04 

2.8382E-06 

6.35 

3.7362E-08 

6.24 

3.7811E-10 

6.62 

3 

B 

1.4253E-04 

6.0409E-06 

4.56 

3.8018E-07 

3.99 

2.3783E-08 

3.99 

A 

1.0433E-04 

4.1471E-07 

7.97 

1.6121E-09 

8.00 

6.2839E-12 

8.00 


U(x,t) = sin(x — t). The numerical errors and order of accuracy before and after post-processing can be 
found in Table 5.4. We clearly observe an accuracy of Ax fc+1 before the post-processing and Ax 2 *+ 2 (one 
order higher than the expected Ax 2fe+1 ) after it, when piecewise P k elements are used. This strongly suggests 
that there is a better negative norm error estimate on the order of Ax 2 * +2 (at least for k = 2 and 3 among 
the tested cases) for piecewise P k elements to this problem. We have used multiple precision FORTRAN 
to compute part of this problem, with an effective precision higher than the standard double precision, to 
overcome some of the round-off problems for this problem. If such multiple precision procedure is not used, 
then the post-processed error can only go down to around 10 -8 to 10 -9 level and cannot decrease further. 
The result in Table 5.4 clearly verifies both the stability and accuracy of the local discontinuous Galerkin 
method developed in section 4 and also the effectiveness of the post-processing technique for such methods 
applied to this type of equations. 


6. Concluding remarks. We have reviewed local discontinuous Galerkin methods for solving PDEs 
with first, second and third spatial derivatives and developed new local discontinuous Galerkin methods 
for solving PDEs with fourth derivatives (the time dependent bi-harmonic equations) and those with fifth 
derivatives. We have designed the numerical fluxes and proved that such methods are L 2 stable and satisfy 
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cell entropy inequalities for general nonlinear PDEs. Although the discussion is concentrate in ID because 
of space limitation, the results can be easily generalized to multiple spatial dimensions. Error estimates 
for linear equations similar to those offered in [11] and [20] can be obtained along similar lines but are not 
given in this paper. PDEs with even higher spatial derivatives can also be handled along similar lines. We 
also apply a post processing technique and demonstrate that the accuracy can be improved from Ax 
to Ax 2k+l or higher when piecewise polynomials of degree k are used. In future work we will investigate 
efficient time marching of these local discontinuous Galerkin methods without sacrificing its local property 
and parallel flexibility, and perform more numerical experiments with physically interesting problems. 
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